General Disclaimer 


One or more of the Following Statements may affect this Document 


• This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 


• This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 


• This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 


• This document is paginated as submitted by the original source. 


• Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 


Produced by the NASA Center for Aerospace Information (CASI) 



NASA TM X-52653 


NASA TECHNICAL 

MEMORANDUM 


rr\ 

LT\ 

SO 

OsJ 

LP\ 

I 

X 




/ViO/ 

CONFORMAL MAPPING PROCEDURE FOR TRANSIENT AND 
STEADY- STATE TWO-DIMENSIONAL SOLIDIFICATION 


by R, Siegel, M, E. Goldstein, and J. M Savino 
Lewis Research Center 
Cleveland, Ohio 


TECHNICAL PAPER proposed for presentation at 
Fourth International Heat Transfer Conference 
Versailles/Paris, August 31 -September 5, 1970 


il7QL - ...2 -3 284 

_Z2 


(TMRUI 


/ 


t IPACUI 

' (NASA CR OR T MX OR AO NUMBLH) 


<33 


ICAVCCtORYI 


oc 


CONFORMAL MAPPING PROCEDURES FOR TRANSIENT AND STEADY-STATE 


TWO-DIMENSIONAL SOLIDIFICATION 

by R Siegel, M. E. Goldstein, and J. M. Savmo 

Lewis Research Center 
Cleveland, Ohio 


TECHNICAL PAPER proposed for presentation at 

Fourth International Heat Transfer Conference 
Versailles/Paris, France, August 31 -September 5, 1970 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


197 


CONFORMAL MAPPING PROCEDURE FOR TRANSIENT AND STEADY -STATE 
TWO-DIMENSIONAL SOLIDIFICATION 

R. Siegel, M. E. Golastein, and J. M. Savino 
NAJ3A Lewis Research Center, Cleveland, Ohio 


Abstract 

A conformal mapping method was developed for analyzing two-dimensional tran- 
sient and steady-state solidification problems. The method was applied to the 
solidification which takes place on a cold plate of finite width immersed in a 
flowing liquid and to the solidification inside of a cooled rectangular channel 
which contains a warm flowing liquid. The transient and steady-state shapes of 
the frozen regions are investigated. 


I NTRODUCTION 

A method was developed for solving two-dimensional transient and steady state 
solidification problems. The method is applicable to the case where a warm 
liquid at a temperature above its equilibrium freezing point flows steadily 
over a surface which is cooled below the freezing point. This may occur, for 
example, inside the conduits of certain rectangular heat exchangers. The 
method is applied to two specific cases which are illustrated in figures 1 
and 2. The first of these consists of the frozen region formed on a cooled 
plate immersed in a warm flowing liquid. The second consists of the frozen 
region which forms inside of a rectangular channel when the channel walls 
are maintained at a constant temperature which is below the freezing tem- 
perature of the liquid. 

In general, the flowing liquid supplies energy by convection to the solid-liquid 
interface. The shape of the frozen region adjusts so that this energy along 
with the latent heat of fusion, which is released in the transient situation, 
can be removed by conduction through the frozen region to the cold boundaries. 

In the transient situation there is also internal energy removed as the frozen 
material is cooled below its freezing point. In the present analysis this en- 
ergy of subcooling is neglected. This assumption is a reasonable one to make 
because in a great many solidification problems, the latent heat released at 
the solid-liquid interface is much greater than the energy of subcooling. 

It is also assumed that the solid-liquid interface is at the equilibrium 
freezing temperature. 


SYMBOLS 

A dimensionless half v^dth, (ha/k)[(tj -tf)/(t f - t w ) J 
A rj time dependent coefficients in mapping 
a half width of plate; half width of long side of channel 

B dimensionless half width, (bb/k)[(tj -t f )/(t f - t w ) ] 


half width of short side of channel 


intermediate mapping parameters 

complete elliptic integral of the second kind 

normal elliptic integral of the first kind 

heat transfer coefficient from flowing liquid to frozen interface 


r r t + i/y ~ i 

m 'c.y " l)(t + Cg)(t - c^)J 


1/2 


dt 


frozen region in Z-plane 
frozen region in W-plane 

complete elliptic integral of the first kind 
.«/ 2 


do), 


0 , 1 « 2 , • 


“ cos 2n^> 

Vl - 3^ sin 2 o) 

thermal conductivity of solidified materia], 
frozen region in £ plane 
(PAq - A) /in V 1 * 0 2 


outward normal 

heat flow rate through frozen layer per unit length 

position vector to frozen interface 

dimensionless temperature, (t - t w )/(t f - t w ) 

temperature 

freezing temperature 

liquid temperature 

surface temperature of cold plate or channel wall 
analytic function, -T + i¥ 
dimensionless coordinates, (x/a)A, (y/a)A 
Cartesian coordinates in physical plane 
dimensionless complex physical plane, X + iY 
physical plane, x + iy 

time dependent coefficients in mapping equation 
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3 time dependent parameter in mapping 

^initial initial value of 3 

frozen region in ft -plane 
n 

r n defined by ^ r k r n-k “ a n* n - 0, 1, 2, . . . 

k=0 

I Kronecker delta 

quantity defined as [-(dT/dx) + i^T/ctt)] -1 
dimensionless time (h‘ /kPA) - t r )‘ /(tj. - t w )]tf 
0 time 

¥ imaginary part of W 

A latent heat of fusion 

| parametric variable 

p density of solidified material 

ft intermediate mapping plane 

a) argument in ft -plane 

Subscript: 

s on frozen interface 

Superscript: 

ss steady state 


GENERAL ANALYSIS 


According to the model adopted the solid-liquid interface is at a constant 
(with both time and position) temperature t^. Since the shape of this 
interface is unknown it is necessary to specify an additional boundary con- 
dition along it. Assume that the heat transfer coefficient h at the solid 
liquid interface is constant. Then at steady state the heat flux into the 
frozen region kn*Vt is uniform along the interface and equal to the convec- 
tive heat supplied by the flowing liquid h(tj - t~) . During the transient, 
however, the rate of freezing is in general nonuniform along the interface 
and the heat flux entering the frozen material is an unknown function of 
position and time which is determined from the condition 

kn*Vt - h(tj - t f ) = pAn (l) 
that results from applying an energy balance at the solid-liquid interface. 

It is convenient to introduce a nondimens ional temperature T defined by 
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T » (t - - t w ). All lengths are nondimetsionalized by 

k./h(tj» - ty)/(tj - ) and the time i3 made nondinensional by 

(kpA/h^) [(t.- - t w )/(tj - tf)j* The dimensional quantities are denoted by 
lower case letters and the dimensionless quantities are denoted by the corre- 
sponding capital letters. 

With the subcooling neglected the heat flow in the solidified region is 
governed by the two-dimensional Laplace's equation (in normalized coordinates), 
that is, at each instant of time the temperature T within the solidified re- 
gion is a harmonic function of position. Hence, let K be the harmonic func- 
tion which is conjugate to -T. Then the complex function W * -T + i^ is a 
function of time and at each fixed instant of time is an analytic function of 
the complex variable Z = X + iY. In view of this we use the notation cW/r)Z 
to denote the ordinary derivative of the analytic function V/ with respect 
to the complex variable Z. 


It is convenient to introduce the complex variable £ defined by 

„ dZ 
S = 3iJ 


( 2 ) 


Clearly, at each instant of time £ is an analytic function of the complex 
variable Z. The function ^ is related to the reciprocal of the complex 
temperature gradient in the frozen region since 


1 dw bT . dT 

l = 3z = ‘ 3x + 1 57 


(3) 


At each instant of time the functions W and £ may be thought of as map- 
pings of the instantaneous frozen region 3© in the physical plane into a 
region J© in the complex W-plane and a region I© in the complex t -plane, 
respectively. Specifying boundary conditions on the functions W and £ 
along the complete boundary of I© is equivalent to specifying the shapes 
of the regions J© and I©. Once the shapes of J© and I© are known, 
it is possible, at least in principle, to introduce an intermediate variable 
ft, choose a certain region P in the ft-plane and then to find the functions 
ft -» W and ft -*■ £ which map T into J© and I©, respectively. When these 
functions are known the integral 

/ ^W / V 

(; ^ dft + function of time, ( 4 ) 

obtained by integrating equation (2), can be evaluated. The integral and the 
known function ft -*• W relate V/ to Z through the parametric variable ft. 
Hence, at each instant of time, the temperature is known at each point of the 
region I© in the physical plane. Since the solid-liquid interface corre- 
sponds to a particular isotherm this correspondence determines the shape of 
the frozen region in the physical plane. Thus, on^e the shapes of the re- 
gions J© and I© are known the solution to the problem can be found. 


Some important differences between the steady state and transient cases 
should be noted. For the steady state cases the regions in the £ and W 
planes can be determined directly from the boundary conditions. This is be- 
cause the uniformity of the heat flux at the solid-liquid interface implies 
that |£| is constant there. In the transient freezing problem the shapes 
of the regions J© and I© change with time, however, it is convenient to 
fix T so that its shape and size are independent of time. Also, part of the 
boundary of I© is unknown and must be determined by solving a nonlinear 
equation. If the region P is suitably chosen; however, the mapping ft £ 
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can be represented by a Taylor series with real coefficients which are func- 
tions of time only. These coefficients can be determined by substituting this 
series into the boundary condition ( 2 ). 


ANALYSIS FCR FREEZING ON PLATE 


The method is best Illustrated by considering the situation depicted in fig- 
ure 1. The cross section of the frozen layer configuration is shown in non- 
dimensional coordinates in figure 3. The nondimens ional boundary conditions 
are all shown on the figure except for the one given by equation (l). Th- 
boundary conditions expressed in terras of the complex variables W and f are 


f e W(Z,0) 
ft <’ W(Z,0) 


-1; Z C FAB 
0; ZC Sc 


JmU(Z ,0) 
Z,0) 


= function of time; 
= function of time; 



Jm^( Z,B) 


ft 4> t (Z,G) 


0 

0 


Z C CB 

7,c? 


Z c EDO 


(3) 

( 6 ) 

(7) 

(e) 


gg 5 * ?(z,©) - U(z,0)| zcfab (9) 

Equations (5) and ( 6 ) show that at each instant of time the region Iq in 
the physical plane maps into the rectangular region Jq as indicated in fig- 
ure 4. The height of the rectangle Jq varies with time and must be deter- 
mined from the soiution to the problem. Equations (7) and ( 8 ) show that the 
region in the physical plane maps into the region 1 ^ in the t; -plane 

as sh-*wn in figure 5. In the transient case the shape of the curve BAF is 
unknown and must be determined by applying equation (9). However, for steady 
state equation (9) becomes |^| =1 and this shows that liA^* is then a semi- 
circle. 


The region T in the ft-plane is chosen to be the semicircular region of unit 
radius shown in figure 6 . Notice that at steady state the mapping Q -* £ 
merely involves multiplication by a negative constant. In the transient situa- 
tion, however, the mapping function is unknown, but an examination of figures 
5 and 6 shows that we may expect this mapping to be continuous on the bounda- 
ries of F since these are no singularities which can occur there. Hence, 
the mapping function can be expanded in a Taylor series about the origin 
which can be expected to converge on the boundaries of T. It also follows 
from figures 5 and 6 that this series must have the form 

C(L,0) = -k(^/i - S 2> ) o n (2 2n+1 (10) 

where the unknown functions of time p and a Q for n = 0 , 1 , 2 , . . . 
are real valued. Elementary mapping techniques can be appD led to show that 
the function which maps T into J 0 in the manner indicated in figures 4 
and 6 is defined by: 
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AW 


t (V~- ") yitl ) - (i - r ; )0 - 


(u) 


Substituting equations (10) and (11 ) into equation (4) and choosing the ori- 
gin oT the coordinate system in the physical plane to be at point D yields 
after ner forming the Integration 


PA 0 -A 

Z( ft,0) * In 

lnVl-P 2 


3 yk(ft) + (l + n 2 ) + (l -3 2 )(i - 2 ) 


vr- 



+ VX(»:) 

n» 


2n 


( 12 ) 

2 2 

with X(n) u (l + i r) - (l - P' )(l - i!^) and the ore related to the 

°n ^ 


«*» “2 ( -) Cl ‘ ( ' 1)J(1 • p2) ] ( n + 1 - 2 ) Vl-j +6 nQp| 
>0 


PA 0 - A 


Ln Vl“" 


, n * 0,1,2, . . . 
(13) 


To determine 3 and the (and hence, the a^) as functions of time, equa- 
tions (13), (12), and (10) are inserted into the boundary condition (9). 

Since ft = e iG (0 < u> < jt) fc: ft C ® the resulting expression involves 

sines and cosines of a). To eliminate the o> dependence of this expression 
it i nultiplied through by cos(2po>) for p » 0, 1, 2, . . . (the restric- 
tion this subset of the complete set of sines and cosines is diet ;ed by 
syiametiy requirements) and integrated oetween (u = 0 and it/2. UpoL per- 
forming these operations we obtain the following infinite set of firs- order 
ordinary differential equations which determine 3 and the A n as functions 
of time. 




and 
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j(r) 

J n,k,p 


rain[2,r) 


E 


2-6 


rl 


(-D J [(-D rJ - (1 * 0 ">M£!.i.k-l),p + b nO»4° 


( 0 ) 

P 


r =1,2,7 and n, k, p ■ 0 , 1, 2, 


• • • 


The Initial values of 3 and A n for these differential equations are cho- 
sen so that the desired initial configuration of the frozen region is given 
by equation (12). 


It is not hard to show that the heat flow through the frozen region can be 
computed from 


2k(t f ‘- V 


kVi - P 2 


(15) 


nee 3 is known. The shape of the frozen region can be computed from 

equal Ion (12) vltb • A 


RESULTS FOR STEADY STATE 


The results for steady state can be obtained as a special case of the pre- 
ceding analysis by letting the time dependent coefficients be zero and 

letting 3 be independent of time. In this case the shape of the solid- 
liquid interface of the frozen layer is given parametrically by 


y 38 

s 


Y ss = 
s 


In 


In 


VTT 


Vi - 


In 


p cos o> + Vl - sin 2 ^ 

Vl - V 




0 £ a) £ it/ 2 


r . 

a) + +«rr- L 

- 

-(l - 3 2 )sin a) 


_3 V 

1 - 3 2 sin 2 ^ + cos os 

t 

s. 


• 

> 


> 


-it/2 S tan"' 1 < 0 


(16) 


The boundary condition on the solid-liquid surface now determines 3 in 
terms of the physical quantity A. Thus, 


In \1 - 

Pk(Vi - P 2 ) 


(IV) 


The heat flow tlirough the frozen region is still given by equation (15). 

For a given A as determined by the imposed temperatures and heat transfer 
coefficient, 3 can be found from equation (17). This value of 3 can 
then be used in equation (16) to compute the shape of the frozen region and 
in equation (15) to compute the heat flow through the frozen region. 


An analysis similar to the one discussed above shows that at steady state in 
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the case of freezing inaiae a rectangular duct the ehape of the frozen region 
is given parametrically by 



Tiie constants c and d are given in terras of the physical parameters of the 
problem by 



The heat flow through the frozen region i3 given tv* 



QUASI-STEADY SOLUTIONS 

In most cases a good approximation to the full transient solution discussed 
above can be obtained by setting all the A^s equal to zero and letting 3 
be the only unknown function of time. This approximation amounts to assuming 
that the heat flux is uniform over the solid liquid interface. In this case, 
the boundary condition on the solid-liquid frozen layer surface implies that 
3 is given as a function of time by 



where P.initial * s the value of p at 0=0. This approximation to the 
solution requires that the frozen region pass through a succession of steady 
states during its transient growth. Thus, the initial configuration must be 
choseu to be a particular steady state frozen region shape. The value of 
^initial the va l ue of 3 determined from the steady state analysis which 
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gives the desired initial steady state shape. Setting ^initial ■ 1 corre- 
sponds to starting the transient when there is no frozen layer on the plate. 


RESULTS AND DISCUSSION 


Tn all the full transient solutions which were carried out the initial con- 
figuration of the frozen region was taken to be a steady state configuration, 
hence, the initial conditions on the Ajj and 0 were taken to be ^ « 0 

and 0 = 3< n itial where ^initial determined from the steady state 
analysis in the same way as for the approximate solution discussed above. 

A typical set of transient growth curves for the frozen layer forming on a 
flat plate are shown in figure 7. The results give a qualitative idea of the 
rate at which solidification occurs . The rate is most rapid at early times 
when the frozen region has the least thickness, and hence, ohe least resis- 
tance to heat flowing through it. As the frozen region becomes large com- 
pared with the width of the plate its shape becomes circular tending toward 
the axisyraraetric solution where the heat removal is through a line sink at 
the center of the solidified region. 

Only the steady state case was considered for freezing inside a rectangular 
duct. A typical set of steady state contours of the frozen region are shown 
in figure 8 . The curves on this plot are drawn for various values of the 
controllable physical parameter B = (hb/k)[(t^ - t^J/Ct^. - t w ) ] . For thin 

layers B is large and the layer is of constant thickness except close to 
the corner. As B is decreased (this corresponds to increased cooling) to 
a value of about 2.5 the frozen layer thickness increases fairly uniformly 
around the duct. Then as B is decreased by a very small amount, the thick- 
ness along the short side increases substantially while tnat of the long side 
remains almost constant. For thick frozen regions the interface approaches 
a circular shape. 

One of the interesting aspects of the profiles in figure 8 is that there is 
a minimum value of B equal to about 2. This phenomena occurs for ducts 
of all aspect ratios. For smaller values of B (i.e., larger cooling) there 
are no steady state solutions and the liquid in the duct would freeze com- 
pletely . 

The fact that there is a minimum value of B leads to the most interesting 
feature of figure 8. For each value of B above the minimum there are 
mathematically two possible frozen region configurations. It can be shown, 
however, that the thicker regions are unstable to small disturbances and 
hence, will not occur physically. 




Figure 3. - Dimensionless physical plane. 
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Figure 5. 
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Figure 8. - Steady-state Irozen layer profiles; duct aspect ratio, 
a/b - 3. 


Figure 7. - Transient solidification starting from an initial layer; 
initial “ 0.995. A - 0. 2. 
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